setwd("~/Dropbox/TingleyMildenberger/ContactCongress/Data/BJPSReplication/data")

##### START OF REPLICATION CODE to SHARE ####

#### FUNCTIONS USED IN CODE ####

## Summarizes data compactly for plotting - credit to: http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/#Helper functions
## Gives count, mean, standard deviation, standard error of the mean, and confidence interval (default 95%).
##   data: a data frame.
##   measurevar: the name of a column that contains the variable to be summariezed
##   groupvars: a vector containing names of columns that contain grouping variables
##   na.rm: a boolean that indicates whether to ignore NA's
##   conf.interval: the percent range of the confidence interval (default is 95%)
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  require(plyr)
  
  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  
  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  datac <- ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2(xx[[col]], na.rm=na.rm),
                     mean = mean   (xx[[col]], na.rm=na.rm),
                     sd   = sd     (xx[[col]], na.rm=na.rm)
                   )
                 },
                 measurevar
  )
  
  # Rename the "mean" column    
  datac <- rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  
  # Confidence interval multiplier for standard error
  # Calculate t-statistic for confidence interval: 
  # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)
}

## Extracts legends as  Grob
g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}


# packages etc
library(foreign)
library(psych)
library(plyr)
library(ggplot2)
library(reshape2)
library(sp)
library(rgeos)
library(Hmisc)
library(grid)
library(gridExtra)



#### FIGURE 1 REPLICATION CODE ####

load(file="SSI2014.Rda")

# calculate population means
mean_happening <- mean(results$gw_happening, na.rm=TRUE) 
mean_human <- mean(results$gw_human, na.rm=TRUE) 
mean_scientists <- mean(results$gw_scientists, na.rm=TRUE) 

# convert to factor
results$gw_happening <- as.factor(results$gw_happening)
results$gw_human <- as.factor(results$gw_human)
results$gw_scientists <- as.factor(results$gw_scientists)


## Perception of population opinions versus individual opinion
national_happening <- summarySE(results, measurevar="gwslider1", groupvars=c("gw_happening"), na.rm=TRUE)
national_happening <- national_happening[1:2,]

national_human <- summarySE(results, measurevar="gwslider2", groupvars=c("gw_human"), na.rm=TRUE)
national_human <- national_human[1:2,]

national_scientists <- summarySE(results, measurevar="gwslider3", groupvars=c("gw_scientists"), na.rm=TRUE)
national_scientists <- national_scientists[1:2,]

national_happening$gw_happening <- c("Disagree", "Agree")
national_happening$gw_happening <- as.factor(national_happening$gw_happening)
national_happening$gw_happening = factor(national_happening$gw_happening,levels(national_happening$gw_happening)[c(2,1)])

national_human$gw_human <- c("Disagree", "Agree")
national_human$gw_human <- as.factor(national_human$gw_human)
national_human$gw_human = factor(national_human$gw_human,levels(national_human$gw_human)[c(2,1)])

national_scientists$gw_scientists <- c("Disagree", "Agree")
national_scientists$gw_scientists <- as.factor(national_scientists$gw_scientists)
national_scientists$gw_scientists = factor(national_scientists$gw_scientists,levels(national_scientists$gw_scientists)[c(2,1)])



natHappen <- ggplot(national_happening, aes(x=gw_happening, y=gwslider1)) + 
  geom_abline(intercept=mean_happening*100, slope=0) +
  geom_point()+
  scale_y_continuous(limits=c(40,80)) +
  ylab("") +
  xlab("GW Happening") +
  theme_bw()+
  geom_errorbar(aes(ymin=gwslider1-ci, ymax=gwslider1+ci),
                width=.1,                  
                position=position_dodge(.9))

natHuman <- ggplot(national_human, aes(x=gw_human, y=gwslider2)) + 
  geom_abline(intercept=mean_human*100, slope=0) +
  scale_y_continuous(limits=c(40,80)) +
  geom_point()+
  ylab("") +
  xlab("GW Human-Caused") +
  theme_bw()+
  geom_errorbar(aes(ymin=gwslider2-ci, ymax=gwslider2+ci),
                width=.1,                    # Width of the error bars
                position=position_dodge(.9))

natScientists <- ggplot(national_scientists, aes(x=gw_scientists, y=gwslider3)) + 
  geom_abline(intercept=mean_scientists*100, slope=0) +
  scale_y_continuous(limits=c(40,80)) +
  ylab("") +
  xlab("GW Consensus") +
  geom_point()+
  theme_bw()+
  geom_errorbar(aes(ymin=gwslider3-ci, ymax=gwslider3+ci),
                width=.1,                    # Width of the error bars
                position=position_dodge(.9))

quartz(width=7, height=5.5)

grid.arrange(arrangeGrob(natHappen, natHuman, natScientists, ncol=2),
             left=textGrob("Estimated fraction of US population who agree", rot = 90, vjust = 1),
             bottom=textGrob("Respondent's Personal Belief",just="bottom"))

#dev.copy2pdf(file="../ContactCongress/Paper/Figures/perceptionsbybeliefSSI2014.pdf")


#### FIGURE 2 ####

load(file="ssichina2015.Rda")

# find population means for reference
mean_chinahappening2 <- mean(as.numeric(as.character(china.results$gw_happening2)), na.rm=TRUE) 
mean_chinatreaty2 <-  mean(as.numeric(as.character(china.results$gw_treaty2)), na.rm=TRUE) 
mean_chinascientists2 <-  mean(as.numeric(as.character(china.results$gw_scientists2)), na.rm=TRUE) 
mean_chinaprice2 <-  mean(as.numeric(as.character(china.results$gw_price2)), na.rm=TRUE) 

# convert to factor
china.results$gw_happening2 <- as.factor(china.results$gw_happening2)
china.results$gw_treaty2 <- as.factor(china.results$gw_treaty2)
china.results$gw_scientists2 <- as.factor(china.results$gw_scientists2)
china.results$gw_price2 <- as.factor(china.results$gw_price2)

china_happening_chinaaboutchina2 <- summarySE(china.results[which(china.results$v116=="cc1|cc2a|cc2b|cc3"),], measurevar="cc2a_1_1_TEXT", groupvars=c("gw_happening2"), na.rm=TRUE)
china_happening_chinaaboutchina2 <- china_happening_chinaaboutchina2[1:2,]
china_treaty_chinaaboutchina2 <- summarySE(china.results[which(china.results$v116=="cc1|cc2a|cc2b|cc3"),], measurevar="cc2a_2_1_TEXT", groupvars=c("gw_treaty2"), na.rm=TRUE)
china_treaty_chinaaboutchina2 <- china_treaty_chinaaboutchina2[1:2,]
china_scientists_chinaaboutchina2 <- summarySE(china.results[which(china.results$v116=="cc1|cc2a|cc2b|cc3"),], measurevar="cc2a_3_1_TEXT", groupvars=c("gw_scientists2"), na.rm=TRUE)
china_scientists_chinaaboutchina2 <- china_scientists_chinaaboutchina2[1:2,]
china_price_chinaaboutchina2 <- summarySE(china.results[which(china.results$v116=="cc1|cc2a|cc2b|cc3"),], measurevar="cc2a_4_1_TEXT", groupvars=c("gw_price2"), na.rm=TRUE)
china_price_chinaaboutchina2 <- china_price_chinaaboutchina2[1:2,]

### graph results
china_happening_chinaaboutchina2$gw_happening <- c("Disagree", "Agree")
china_happening_chinaaboutchina2$gw_happening <- as.factor(china_happening_chinaaboutchina2$gw_happening)
china_happening_chinaaboutchina2$gw_happening = factor(china_happening_chinaaboutchina2$gw_happening,levels(china_happening_chinaaboutchina2$gw_happening)[c(2,1)])

china_treaty_chinaaboutchina2$gw_treaty <- c("Disagree", "Agree")
china_treaty_chinaaboutchina2$gw_treaty <- as.factor(china_treaty_chinaaboutchina2$gw_treaty)
china_treaty_chinaaboutchina2$gw_treaty <- factor(china_treaty_chinaaboutchina2$gw_treaty,levels(china_treaty_chinaaboutchina2$gw_treaty)[c(2,1)])

china_scientists_chinaaboutchina2$gw_scientists <- c("Disagree", "Agree")
china_scientists_chinaaboutchina2$gw_scientists <- as.factor(china_scientists_chinaaboutchina2$gw_scientists)
china_scientists_chinaaboutchina2$gw_scientists <- factor(china_scientists_chinaaboutchina2$gw_scientists,levels(china_scientists_chinaaboutchina2$gw_scientists)[c(2,1)])

china_price_chinaaboutchina2$gw_price <- c("Disagree", "Agree")
china_price_chinaaboutchina2$gw_price <- as.factor(china_price_chinaaboutchina2$gw_price)
china_price_chinaaboutchina2$gw_price <- factor(china_price_chinaaboutchina2$gw_price,levels(china_price_chinaaboutchina2$gw_price)[c(2,1)])


chinamean_happening <- mean(as.numeric(as.character(china.results$gw_happening2)), na.rm=TRUE)
chinamean_treaty <- mean(as.numeric(as.character(china.results$gw_treaty2)), na.rm=TRUE)
chinamean_scientists <- mean(as.numeric(as.character(china.results$gw_scientists2)), na.rm=TRUE)
chinamean_price <- mean(as.numeric(as.character(china.results$gw_price2)), na.rm=TRUE)


chinaHappen2 <- ggplot(china_happening_chinaaboutchina2, aes(x=gw_happening, y=cc2a_1_1_TEXT)) + 
  geom_abline(intercept=chinamean_happening*100, slope=0) +
  geom_point()+
  scale_y_continuous(limits=c(10,100)) +
  ylab("") +
  xlab("GW Happening") +
  theme_bw()+
  geom_errorbar(aes(ymin=cc2a_1_1_TEXT-ci, ymax=cc2a_1_1_TEXT+ci),
                width=.1,                  
                position=position_dodge(.9))

chinaTreaty2 <- ggplot(china_treaty_chinaaboutchina2, aes(x=gw_treaty, y=cc2a_2_1_TEXT)) + 
  geom_abline(intercept=chinamean_treaty*100, slope=0) +
  geom_point()+
  scale_y_continuous(limits=c(10,100)) +
  ylab("") +
  xlab("GW Treaty") +
  theme_bw()+
  geom_errorbar(aes(ymin=cc2a_2_1_TEXT-ci, ymax=cc2a_2_1_TEXT+ci),
                width=.1,                  
                position=position_dodge(.9))

chinaScientists2 <- ggplot(china_scientists_chinaaboutchina2, aes(x=gw_scientists, y=cc2a_3_1_TEXT)) + 
  geom_abline(intercept=chinamean_scientists*100, slope=0) +
  geom_point()+
  scale_y_continuous(limits=c(10,100)) +
  ylab("") +
  xlab("GW Consensus") +
  theme_bw()+
  geom_errorbar(aes(ymin=cc2a_3_1_TEXT-ci, ymax=cc2a_3_1_TEXT+ci),
                width=.1,                  
                position=position_dodge(.9))

chinaPrice2 <- ggplot(china_price_chinaaboutchina2, aes(x=gw_price, y=cc2a_4_1_TEXT)) + 
  geom_abline(intercept=chinamean_price*100, slope=0) +
  geom_point()+
  scale_y_continuous(limits=c(10,100)) +
  ylab("") +
  xlab("GW Price") +
  theme_bw()+
  geom_errorbar(aes(ymin=cc2a_4_1_TEXT-ci, ymax=cc2a_4_1_TEXT+ci),
                width=.1,                  
                position=position_dodge(.9))

quartz(width=7, height=5.5)

grid.arrange(arrangeGrob(chinaHappen2, chinaTreaty2, chinaScientists2, chinaPrice2, ncol=2),
             left=textGrob("Estimated fraction of the Chinese population who agree", rot = 90, vjust = 1),
             bottom=textGrob("Respondent's Personal Belief",just="bottom"))

#dev.copy2pdf(file="../ContactCongress/Paper/Figures/chinageographysplits.pdf")

#### FIGURE 3 ####
load(file="MTurkMay2014.Rda")

geoHappening <- summarySE(exp_May, measurevar="geo_happening", groupvars=c("happenBIN","country"), na.rm=TRUE)
geoHappening2 <- geoHappening[-c(2,3,5,7,8),]
geoHappening3 <- geoHappening[-c(4,5,9),]


geo1 <- ggplot(geoHappening2, aes(x=country, y=geo_happening, shape=happenBIN)) + 
  geom_point(position=position_dodge(.2)) +
  ylab("") +
  xlab("GW Happening") +
  scale_shape_discrete(name="Respondent",
                       breaks=c("0","1"),
                       labels=c("Disagrees","Agrees"), solid=FALSE) +
  theme_bw()+
  scale_y_continuous(limits=c(25,75)) +
  geom_errorbar(aes(ymin=geo_happening-ci, ymax=geo_happening+ci),
                width=.1, position=position_dodge(.2)) +
  theme(legend.position="bottom")

geoTreaty <- summarySE(exp_May, measurevar="geo_treaty", groupvars=c("treatyBIN","country"), na.rm=TRUE)
geoTreaty2 <- geoTreaty[-c(2,3,5,7,8),]
geoTreaty3 <- geoTreaty[-c(4,5,9),]

geo2 <- ggplot(geoTreaty2, aes(x=country, y=geo_treaty, shape=treatyBIN)) + 
  geom_point(position=position_dodge(.2)) +
  ylab("") +
  xlab("GW Treaty") +
  scale_shape_discrete(name="Respondent",
                       breaks=c("0","1"),
                       labels=c("Disagrees","Agrees"), solid=FALSE) +
  theme_bw()+
  scale_y_continuous(limits=c(25,75)) +
  geom_errorbar(aes(ymin=geo_treaty-ci, ymax=geo_treaty+ci),
                width=.1, position=position_dodge(.2)) +
  theme(legend.position="none")

geoConsensus <- summarySE(exp_May, measurevar="geo_consensus", groupvars=c("consensusBIN","country"), na.rm=TRUE)
geoConsensus2 <- geoConsensus[-c(2,3,5,7,8),]
geoConsensus3 <- geoConsensus[-c(4,5,9),]

geo3 <- ggplot(geoConsensus2, aes(x=country, y=geo_consensus, shape=consensusBIN)) + 
  geom_point(position=position_dodge(.2)) +
  ylab("") +
  xlab("GW Consensus") +
  scale_shape_discrete(name="Respondent",
                       breaks=c("0","1"),
                       labels=c("Disagrees","Agrees"), solid=FALSE) +
  theme_bw()+
  scale_y_continuous(limits=c(25,75)) +
  geom_errorbar(aes(ymin=geo_consensus-ci, ymax=geo_consensus+ci),
                width=.1, position=position_dodge(.2)) +
  theme(legend.position="none")

geoPrice <- summarySE(exp_May, measurevar="geo_price", groupvars=c("priceBIN","country"), na.rm=TRUE)
geoPrice2 <- geoPrice[-c(2,3,5,7,8),]
geoPrice3 <- geoPrice[-c(4,5,9),]


geo4 <- ggplot(geoPrice2, aes(x=country, y=geo_price, shape=priceBIN)) + 
  geom_point(position=position_dodge(.2)) +
  ylab("") +
  xlab("GW Price") +
  scale_shape_discrete(name="Respondent",
                       breaks=c("0","1"),
                       labels=c("Disagrees","Agrees"), solid=FALSE) +
  theme_bw()+
  scale_y_continuous(limits=c(25,75)) +
  geom_errorbar(aes(ymin=geo_price-ci, ymax=geo_price+ci),
                width=.1, position=position_dodge(.2)) +
  theme(legend.position="none")

countrylegend <- g_legend(geo1)
geo1 <- geo1 + theme(legend.position="none")

quartz(width=7, height=5.5)
grid.arrange(arrangeGrob(geo1, geo2, geo3, geo4, ncol=2), countrylegend, nrow=2,heights=c(10, 1),
             left=textGrob("Estimated fraction of country's population who agrees", rot=90, vjust=1))
#dev.copy2pdf(file="../ContactCongress/Paper/Figures/chinaUScontrast_USSurvey.pdf")



### FIGURE 4 ####

china_happening_usaboutchina2 <- summarySE(china.results, measurevar="cc2b_1_1_TEXT", groupvars=c("gw_happening2"), na.rm=TRUE)
china_happening_usaboutchina2 <- china_happening_usaboutchina2[1:2,]
china_treaty_usaboutchina2 <- summarySE(china.results, measurevar="cc2b_2_1_TEXT", groupvars=c("gw_treaty2"), na.rm=TRUE)
china_treaty_usaboutchina2 <- china_treaty_usaboutchina2[1:2,]
china_scientists_usaboutchina2 <- summarySE(china.results, measurevar="cc2b_3_1_TEXT", groupvars=c("gw_scientists2"), na.rm=TRUE)
china_scientists_usaboutchina2 <- china_scientists_usaboutchina2[1:2,]
china_price_usaboutchina2 <- summarySE(china.results, measurevar="cc2b_4_1_TEXT", groupvars=c("gw_price2"), na.rm=TRUE)
china_price_usaboutchina2 <- china_price_usaboutchina2[1:2,]

china_happening_usaboutchina2$gw_happening <- c("Disagree","Agree")
china_happening_usaboutchina2$gw_happening <- as.factor(china_happening_usaboutchina2$gw_happening)
china_happening_usaboutchina2$gw_happening = factor(china_happening_usaboutchina2$gw_happening,levels(china_happening_usaboutchina2$gw_happening)[c(2,1)])
china_happening_usaboutchina2$country <- NA
china_happening_usaboutchina2$country <- "US"
names(china_happening_usaboutchina2)[3] <- "geo_happening"

china_happening_chinaaboutchina2$country <- NA
china_happening_chinaaboutchina2$country <- "China"
names(china_happening_chinaaboutchina2)[3] <- "geo_happening"
chinaforeignhappening <- rbind(china_happening_chinaaboutchina2, china_happening_usaboutchina2)

china_treaty_usaboutchina2$gw_treaty <- c("Disagree","Agree")
china_treaty_usaboutchina2$gw_treaty <- as.factor(china_treaty_chinaaboutchina2$gw_treaty)
china_treaty_usaboutchina2$gw_treaty = factor(china_treaty_usaboutchina2$gw_treaty,levels(china_treaty_usaboutchina2$gw_treaty)[c(2,1)])
china_treaty_usaboutchina2$country <- NA
china_treaty_usaboutchina2$country <- "US"
names(china_treaty_usaboutchina2)[3] <- "geo_treaty"

china_treaty_chinaaboutchina2$country <- NA
china_treaty_chinaaboutchina2$country <- "China"
names(china_treaty_chinaaboutchina2)[3] <- "geo_treaty"
chinaforeigntreaty <- rbind(china_treaty_chinaaboutchina2, china_treaty_usaboutchina2)


china_scientists_usaboutchina2$gw_scientists <- c("Disagree","Agree")
china_scientists_usaboutchina2$gw_scientists <- as.factor(china_scientists_usaboutchina2$gw_scientists)
china_scientists_usaboutchina2$gw_scientists = factor(china_scientists_usaboutchina2$gw_scientists,levels(china_scientists_usaboutchina2$gw_scientists)[c(2,1)])
china_scientists_usaboutchina2$country <- NA
china_scientists_usaboutchina2$country <- "US"
names(china_scientists_usaboutchina2)[3] <- "geo_scientists"

china_scientists_chinaaboutchina2$country <- NA
china_scientists_chinaaboutchina2$country <- "China"
names(china_scientists_chinaaboutchina2)[3] <- "geo_scientists"
chinaforeignscientists <- rbind(china_scientists_chinaaboutchina2, china_scientists_usaboutchina2)

china_price_usaboutchina2$gw_price <- c("Disagree","Agree")
china_price_usaboutchina2$gw_price <- as.factor(china_price_usaboutchina2$gw_price)
china_price_usaboutchina2$gw_price  = factor(china_price_usaboutchina2$gw_price ,levels(china_price_usaboutchina2$gw_price )[c(2,1)])
china_price_usaboutchina2$country <- NA
china_price_usaboutchina2$country <- "US"
names(china_price_usaboutchina2)[3] <- "geo_price"

china_price_chinaaboutchina2$country <- NA
china_price_chinaaboutchina2$country <- "China"
names(china_price_chinaaboutchina2)[3] <- "geo_price"
chinaforeignprice <- rbind(china_price_chinaaboutchina2, china_price_usaboutchina2)


chinaUShappen <- ggplot(chinaforeignhappening, aes(x=country, y=geo_happening, shape=gw_happening)) + 
  geom_point(position=position_dodge(.2)) +
  ylab("") +
  xlab("GW Happening") +
  scale_shape_discrete(name="Respondent",
                       breaks=c("0","1"),
                       labels=c("Disagrees","Agrees"), solid=FALSE) +
  theme_bw()+
  scale_y_continuous(limits=c(10,75)) +
  geom_errorbar(aes(ymin=geo_happening-ci, ymax=geo_happening+ci),
                width=.1, position=position_dodge(.2)) +
  theme(legend.position="bottom")

chinaUStreaty <- ggplot(chinaforeigntreaty, aes(x=country, y=geo_treaty, shape=gw_treaty)) + 
  geom_point(position=position_dodge(.2)) +
  ylab("") +
  xlab("GW Treaty") +
  scale_shape_discrete(name="Respondent",
                       breaks=c("0","1"),
                       labels=c("Disagrees","Agrees"), solid=FALSE) +
  theme_bw()+
  scale_y_continuous(limits=c(10,75)) +
  geom_errorbar(aes(ymin=geo_treaty-ci, ymax=geo_treaty+ci),
                width=.1, position=position_dodge(.2)) +
  theme(legend.position="bottom")


chinaUSscientists <- ggplot(chinaforeignscientists, aes(x=country, y=geo_scientists, shape=gw_scientists)) + 
  geom_point(position=position_dodge(.2)) +
  ylab("") +
  xlab("GW Scientists") +
  scale_shape_discrete(name="Respondent",
                       breaks=c("0","1"),
                       labels=c("Disagrees","Agrees"), solid=FALSE) +
  theme_bw()+
  scale_y_continuous(limits=c(10,75)) +
  geom_errorbar(aes(ymin=geo_scientists-ci, ymax=geo_scientists+ci),
                width=.1, position=position_dodge(.2)) +
  theme(legend.position="bottom")

chinaUSprice <- ggplot(chinaforeignprice, aes(x=country, y=geo_price, shape=gw_price)) + 
  geom_point(position=position_dodge(.2)) +
  ylab("") +
  xlab("GW Price") +
  scale_shape_discrete(name="Respondent",
                       breaks=c("0","1"),
                       labels=c("Disagrees","Agrees"), solid=FALSE) +
  theme_bw()+
  scale_y_continuous(limits=c(10,75)) +
  geom_errorbar(aes(ymin=geo_price-ci, ymax=geo_price+ci),
                width=.1, position=position_dodge(.2)) +
  theme(legend.position="bottom")

quartz(width=7, height=5.5)

grid.arrange(arrangeGrob(chinaUShappen, chinaUStreaty, chinaUSscientists, chinaUSprice, ncol=2),countrylegend,nrow=2,heights=c(10, 1),
             left=textGrob("Estimated fraction of country's population who agrees", rot = 90, vjust = 1))
             
# dev.copy2pdf(file="../ContactCongress/Paper/Figures/chinaUScontrast_ChinaSurvey.pdf")

#### FIGURE 5 ####

# seperate script 

#### FIGURE 6 ####

# See seperate script

#### FIGURE 7 ####

load(file="legislativestaffer_aug2016.Rda")

co2_national <- summarySE(data=staffer.extract, measurevar=c("nationalrep_co2"), groupvars=c("personal_co2"), na.rm=TRUE)
co2_national <- co2_national[-3,]
co2_national$personal_co2 <- as.factor(co2_national$personal_co2)
co2_national$personal_co2 <- revalue(co2_national$personal_co2, c("0"="Disagree", "1"="Agree"))

quartz(width=7, height=5)

natCO2_staffer <- ggplot(co2_national, aes(x=personal_co2, y=nationalrep_co2)) + 
  geom_point()+
  scale_y_continuous(limits=c(25,90)) +
  ggtitle("Regulate carbon dioxide")+
  geom_hline(yintercept = 63, linetype="dashed") +
  ylab("Staffer Est. National Agreement") +
  xlab("Staffer Personal Belief") +
  theme_bw()+
  geom_errorbar(aes(ymin=nationalrep_co2-ci, ymax=nationalrep_co2+ci),
                width=.1,                  
                position=position_dodge(.9))
natCO2_staffer
#dev.copy2pdf(file="PaperFigures/stafferegocentricbias.pdf")

#### FIGURE 8 ####
 
# See seperate script


#### ONLINE APPENDIX: FIGURE 1 ####

#sort sample into high and low education to see if there are differences in egocentric bias

lowedu <- results[which(results$educ==1|results$educ==2),]
mededu <- results[which(results$educ==3|results$educ==4|results$educ==5),]
highedu <- results[which(results$educ==6|results$educ==7),]

lowedu_happening <- summarySE(lowedu, measurevar="gwslider1", groupvars=c("gw_happening"), na.rm=TRUE)
lowedu_happening <- lowedu_happening[1:2,]

mededu_happening <- summarySE(mededu, measurevar="gwslider1", groupvars=c("gw_happening"), na.rm=TRUE)
mededu_happening <- mededu_happening[1:2,]

highedu_happening <- summarySE(highedu, measurevar="gwslider1", groupvars=c("gw_happening"), na.rm=TRUE)
highedu_happening <- highedu_happening[1:2,]

lowHappen <- ggplot(lowedu_happening, aes(x=gw_happening, y=gwslider1)) + 
  geom_abline(intercept=mean_happening*100, slope=0) +
  geom_point()+
  scale_y_continuous(limits=c(40,80)) +
  ylab("") +
  xlab("Low education") +
  theme_bw()+
  geom_errorbar(aes(ymin=gwslider1-ci, ymax=gwslider1+ci),
                width=.1,                  
                position=position_dodge(.9))

medHappen <- ggplot(mededu_happening, aes(x=gw_happening, y=gwslider1)) + 
  geom_abline(intercept=mean_happening*100, slope=0) +
  geom_point()+
  scale_y_continuous(limits=c(40,80)) +
  ylab("") +
  xlab("Medium education") +
  theme_bw()+
  geom_errorbar(aes(ymin=gwslider1-ci, ymax=gwslider1+ci),
                width=.1,                  
                position=position_dodge(.9))

highHappen <- ggplot(highedu_happening, aes(x=gw_happening, y=gwslider1)) + 
  geom_abline(intercept=mean_happening*100, slope=0) +
  geom_point()+
  scale_y_continuous(limits=c(40,80)) +
  ylab("") +
  xlab("High education") +
  theme_bw()+
  geom_errorbar(aes(ymin=gwslider1-ci, ymax=gwslider1+ci),
                width=.1,                  
                position=position_dodge(.9))

grid.arrange(arrangeGrob(lowHappen, medHappen, highHappen, ncol=2),
             left=textGrob("Estimated fraction of US population who agree", rot = 90, vjust = 1),
             bottom=textGrob("Respondent's Personal Belief",just="bottom"))

####  ONLINE APPENDIX: FIGURES 2 & 3 ####
### Now to look at differences across partisan answers

poliHappening <- summarySE(exp_May, measurevar="poli_happening", groupvars=c("partyhappening","partyassign"), na.rm=TRUE)
poliHappening_Rep <- poliHappening[c(2,5,7),]
poliHappening_Dem <- poliHappening[c(1,4,6),]

poli1_Rep <- ggplot(poliHappening_Rep, aes(x=partyassign, y=poli_happening, shape=partyhappening)) + 
  geom_point(position=position_dodge(.2)) +
  scale_y_continuous(limits=c(15, 60)) +
  ylab("") +
  xlab("GW Happening") +
  scale_x_discrete(labels=c(""))+
  scale_shape_discrete(name="Respondent Party", solid=FALSE)+
  theme_bw()+
  geom_errorbar(aes(ymin=poli_happening-ci, ymax=poli_happening+ci),
                width=.1,                    # Width of the error bars
                position=position_dodge(.2)) +
  theme(legend.position="bottom")

poli1_Dem <- ggplot(poliHappening_Dem, aes(x=partyassign, y=poli_happening, shape=partyhappening)) + 
  geom_point(position=position_dodge(.2)) +
  scale_y_continuous(limits=c(60, 90)) +
  ylab("") +
  xlab("GW Happening") +
  scale_x_discrete(labels=c(""))+
  scale_shape_discrete(name="Respondent Party", solid=FALSE)+
  theme_bw()+
  geom_errorbar(aes(ymin=poli_happening-ci, ymax=poli_happening+ci),
                width=.1,                    # Width of the error bars
                position=position_dodge(.2)) +
  theme(legend.position="none")

poliTreaty <- summarySE(exp_May, measurevar="poli_treaty", groupvars=c("partytreaty","partyassign"), na.rm=TRUE)
poliTreaty_Rep <- poliTreaty[c(2,5,7),]
poliTreaty_Dem <- poliTreaty[c(1,4,6),]


poli2_Rep <- ggplot(poliTreaty_Rep, aes(x=partyassign, y=poli_treaty, shape=partytreaty)) + 
  geom_point(position=position_dodge(.2)) +
  scale_y_continuous(limits=c(15, 60)) +
  ylab("") +
  xlab("GW Treaty") +
  scale_x_discrete(labels=c(""))+
  scale_shape_discrete(name="Respondent Party", solid=FALSE)+
  theme_bw()+
  geom_errorbar(aes(ymin=poli_treaty-ci, ymax=poli_treaty+ci),
                width=.1,                    # Width of the error bars
                position=position_dodge(.2))+
  theme(legend.position="none")

poli2_Dem <- ggplot(poliTreaty_Dem, aes(x=partyassign, y=poli_treaty, shape=partytreaty)) + 
  geom_point(position=position_dodge(.2)) +
  scale_y_continuous(limits=c(60, 90)) +
  ylab("") +
  xlab("GW Treaty") +
  scale_x_discrete(labels=c(""))+
  scale_shape_discrete(name="Respondent Party", solid=FALSE)+
  theme_bw()+
  geom_errorbar(aes(ymin=poli_treaty-ci, ymax=poli_treaty+ci),
                width=.1,                    # Width of the error bars
                position=position_dodge(.2))+
  theme(legend.position="none")

poliConsensus <- summarySE(exp_May, measurevar="poli_consensus", groupvars=c("partyconsensus","partyassign"), na.rm=TRUE)
poliConsensus_Rep <- poliConsensus[c(2,5,7),]
poliConsensus_Dem <- poliConsensus[c(1,4,6),]


poli3_Rep <- ggplot(poliConsensus_Rep, aes(x=partyassign, y=poli_consensus, shape=partyconsensus)) + 
  geom_point(position=position_dodge(.2)) +
  scale_y_continuous(limits=c(15, 60)) +
  ylab("") +
  xlab("GW Consensus") +
  scale_x_discrete(labels=c(""))+
  scale_shape_discrete(name="Respondent Party", solid=FALSE)+
  theme_bw()+
  geom_errorbar(aes(ymin=poli_consensus-ci, ymax=poli_consensus+ci),
                width=.1,                    # Width of the error bars
                position=position_dodge(.2))+
  theme(legend.position="none")

poli3_Dem <- ggplot(poliConsensus_Dem, aes(x=partyassign, y=poli_consensus, shape=partyconsensus)) + 
  geom_point(position=position_dodge(.2)) +
  scale_y_continuous(limits=c(60, 90)) +
  ylab("") +
  xlab("GW Consensus") +
  scale_x_discrete(labels=c(""))+
  scale_shape_discrete(name="Respondent Party", solid=FALSE)+
  theme_bw()+
  geom_errorbar(aes(ymin=poli_consensus-ci, ymax=poli_consensus+ci),
                width=.1,                    # Width of the error bars
                position=position_dodge(.2))+
  theme(legend.position="none")


poliPrice <- summarySE(exp_May, measurevar="poli_price", groupvars=c("partyprice","partyassign"), na.rm=TRUE)
poliPrice_Rep <- poliPrice[c(2,5,7),]
poliPrice_Dem <- poliPrice[c(1,4,6),]

poli4_Rep <- ggplot(poliPrice_Rep, aes(x=partyassign, y=poli_price, shape=partyprice)) + 
  geom_point(position=position_dodge(.2)) +
  scale_y_continuous(limits=c(15, 60)) +
  ylab("") +
  xlab("GW Price") +
  scale_x_discrete(labels=c(""))+
  scale_shape_discrete(name="Respondent Party", solid=FALSE)+
  theme_bw()+
  geom_errorbar(aes(ymin=poli_price-ci, ymax=poli_price+ci),
                width=.1,                    # Width of the error bars
                position=position_dodge(.2))+
  theme(legend.position="none")

poli4_Dem <- ggplot(poliPrice_Dem, aes(x=partyassign, y=poli_price, shape=partyprice)) + 
  geom_point(position=position_dodge(.2)) +
  scale_y_continuous(limits=c(60, 90)) +
  ylab("") +
  xlab("GW Price") +
  scale_x_discrete(labels=c(""))+
  scale_shape_discrete(name="Respondent Party", solid=FALSE)+
  theme_bw()+
  geom_errorbar(aes(ymin=poli_price-ci, ymax=poli_price+ci),
                width=.1,                    # Width of the error bars
                position=position_dodge(.2))+
  theme(legend.position="none")



partylegend <- g_legend(poli1_Rep) #extract from any of the individual figures
poli1_Rep <- poli1_Rep + theme(legend.position="none")

quartz(width=7, height=5.5)

grid.arrange(arrangeGrob(poli1_Dem, poli2_Dem, poli3_Dem, poli4_Dem, ncol=2), partylegend, nrow=2,heights=c(10, 1),
             left=textGrob("Fraction of Democrats who agree (est.)", rot=90, vjust=1))
#dev.copy2pdf(file="../ContactCongress/Paper/Figures/democraticsplits.pdf")

grid.arrange(arrangeGrob(poli1_Rep, poli2_Rep, poli3_Rep, poli4_Rep, ncol=2), partylegend, nrow=2,heights=c(10, 1),
             left=textGrob("Fraction of Republicans who agree (est.)", rot=90, vjust=1))

#dev.copy2pdf(file="../ContactCongress/Paper/Figures/republicansplits.pdf")


#### ONLINE APPENDIX: FIGURE 4 ####

geo31 <- ggplot(geoHappening3, aes(x=country, y=geo_happening, shape=happenBIN)) + 
  geom_point(position=position_dodge(.2)) +
  ylab("") +
  xlab("GW Happening") +
  scale_shape_discrete(name="Respondent",
                       breaks=c("0","1"),
                       labels=c("Disagrees","Agrees"), solid=FALSE) +
  theme_bw()+
  scale_y_continuous(limits=c(25,85)) +
  geom_errorbar(aes(ymin=geo_happening-ci, ymax=geo_happening+ci),
                width=.1, position=position_dodge(.2)) +
  theme(legend.position="bottom")



geo32 <- ggplot(geoTreaty3, aes(x=country, y=geo_treaty, shape=treatyBIN)) + 
  geom_point(position=position_dodge(.2)) +
  ylab("") +
  xlab("GW Treaty") +
  scale_shape_discrete(name="Respondent",
                       breaks=c("0","1"),
                       labels=c("Disagrees","Agrees"), solid=FALSE) +
  theme_bw()+
  scale_y_continuous(limits=c(25,85)) +
  geom_errorbar(aes(ymin=geo_treaty-ci, ymax=geo_treaty+ci),
                width=.1, position=position_dodge(.2)) +
  theme(legend.position="none")



geo33 <-ggplot(geoConsensus3, aes(x=country, y=geo_consensus, shape=consensusBIN)) + 
  geom_point(position=position_dodge(.2)) +
  ylab("") +
  xlab("GW Consensus") +
  scale_shape_discrete(name="Respondent",
                       breaks=c("0","1"),
                       labels=c("Disagrees","Agrees"), solid=FALSE) +
  theme_bw()+
  scale_y_continuous(limits=c(25,85)) +
  geom_errorbar(aes(ymin=geo_consensus-ci, ymax=geo_consensus+ci),
                width=.1, position=position_dodge(.2)) +
  theme(legend.position="none")

geo34 <- ggplot(geoPrice3, aes(x=country, y=geo_price, shape=priceBIN)) + 
  geom_point(position=position_dodge(.2)) +
  ylab("") +
  xlab("GW Price") +
  scale_shape_discrete(name="Respondent",
                       breaks=c("0","1"),
                       labels=c("Disagrees","Agrees"), solid=FALSE) +
  theme_bw()+
  scale_y_continuous(limits=c(25,85)) +
  geom_errorbar(aes(ymin=geo_price-ci, ymax=geo_price+ci),
                width=.1, position=position_dodge(.2)) +
  theme(legend.position="none")

country3legend <- g_legend(geo31)
geo31 <- geo31 + theme(legend.position="none")

quartz(width=7, height=5.5)
grid.arrange(arrangeGrob(geo31, geo32, geo33, geo34, ncol=2), country3legend, nrow=2,heights=c(10, 1),
             left=textGrob("Estimated fraction of country's population who agrees", rot=90, vjust=1))
#dev.copy2pdf(file="../ContactCongress/Paper/Figures/extendedcontrast_USSurvey.pdf")




